%!TEX root = all.tex

\section{Experiments}
As already mentioned before, we decided to verify the performance of the ICE algorithm
implemented by the C-SAFE (\url{http://www.csafe.utah.edu}) research group at the University
of Utah. The multimaterial ICE algorithm is utilized to simulate explosions, fires and other fluid phenomena.
It is a cell-centered, finite volume algorithm developed and described by
Kashiwa, et. al. \cite{R9}. It uses a gradient limiter \cite{R12} to suppress unnatural oscillations
introduced by the higher-order numerical methods. The effect of the gradient limiter on
the order of accuracy will be discussed below.
\\
\\
The ICE code contains several modules that are designed to solve various sub-problems.  We decided to focus our initial
research on the Advect and Advance in Time (AAT) module.  This selection was undertaken for the following reasons:

\begin{enumerate}
\item The AAT module is a crucial part of ICE simulations; it is invoked multiple times
during full scale fluid simulations.  Therefore, its numerical accuracy is essential.
\item The advection operator is defined by one governing equation, which makes the
multivariate MGS implementations memory efficient on various grid sizes.
\end{enumerate}

\noindent To isolate the Advect and Advance in Time module, the transport of a passive
scalar was employed as a verification experiment. The initial condition of the passive scalar were
defined using a bell-shaped exponential function. The 3D
governing equation for the experiment was as follows:

\begin{equation}
\frac{d}{dt}(PS(x, y, z, t)) + U_x\frac{d}{dx}(PS(x, y, z, t)) + U_y\frac{d}{dy}(PS(x, y, z, t)) + U_z\frac{d}{dz}(PS(x, y, z, t)) = g(x, y, z, t)
\end{equation}
where,
\begin{align*}
&PS(x,\ y,\ z,\ t) - \text{Passive scalar}\\
&U_x,\ U_y,\ U_z - \text{Constant velocity in x,  y, and z directions correspondingly}\\
&g(x,\ y,\ z,\ t) - \text{Source terms}
\end{align*}

As already mentioned above, the MGS method utilizes approximate discrete solution given by
either experimental data or a computational tool. Since physical
experimental data was not available, we decided to use approximate solution generated by
ICE to build our comparison solution. This completes step 1 of the verification process.
We decided to use natural cubic splines as the basis for our interpolation to construct the exact
solution of equation \ref{eq:maineq2}. Since the differential operator D from equation (1)
corresponding to equation (6) is a degree one function, the analytical solution resulting from 
interpolation of the simulated solution \emph{u} has to be at least \emph{C$^1$} continuous. Since natural
cubic spline interpolation produces \emph{C$^2$} continuous functions, which fits the requirement, it
was chosen as an interpolation technique for computing exact solution \emph{u$_{1}$} and gradients at
the specified points of time and space. These gradients are then used to compute source
terms \emph{g$_{1}$(x,y,z,t)} on the right-hand side of the equation (6). These source terms are then substituted 
for the source terms into ICE. Using these new source terms, equation (2) is solved to generate the solution \emph{u$_{2}$}.
Finally, we compare the computed solution \emph{u$_{2}$} with the generated solution \emph{u$_{1}$} and compute L$_{2}$ and L$_{\infty}$ error, error consistency and order of accuracy.
\\
\\
{\bf Test Parameters}
\begin{enumerate}
\item For our verification experiments we decided to choose the initial distribution of the passive scalar to
be a 3D bellshaped exponential.
\item The computational domain was 1x1x1m$^3$ (from -0.5 to 0.5
meters in each direction). We also chose time step, $\Delta$t = 10$^{-6}$, and the number of time
steps taken was 10.
\item Tests were conducted with resolutions 40x40x40, 80x80x80 and 160x160x160 which correspond to $\Delta$x, $\Delta$y and $\Delta$z equal to 1x10$^{-2}$, 0.5x10$^{-2}$ and 0.25x10$^{-2}$ respectively.
\item The source terms \emph{g$_{1}$(x,y,z,t)} are computed using the MGS method from two different input sources:
\begin{compactenum}
\item Exact analytical solution;
\item ICE generated solution.
\end{compactenum}
\item ICE allows the user to specify the order of the discretization scheme used.  We decided to verify its performance when we select both first and second orders.
\item As already mentioned before, ICE utilizes a gradient limiter to suppress non-physical
oscillations in high gradient regions. The limiter
implemented in ICE is computed using a Van Leer type method \cite{R9}\cite{R11}:
\begin{equation*}
\alpha_{j} = min(1,\alpha_{j_{min}},\alpha_{j_{max}})
\end{equation*}
\begin{align*}
&\alpha_{j_{min}} = max\left(0,\frac{u_{min} - u_{j}}{min[u_{v}] - u_{j}}\right), \alpha_{j_{max}} = max\left(0, \frac{u_{max} - u_{j}}{max[u_{v}] - u_{j}}\right),\\
\end{align*}
where,
\begin{align*}
&u_{j} - \text{value of the solution at the cell center;}\\
&u_{v} - \text{value of the solution at the cell vertices;}\\
&u_{min},u_{max} - \text{min and max values of the solution at the surrounding cell centers.}\\
\end{align*}
The values of the gradient limiters are used to bound the values of gradients:
\begin{equation*}
(\Delta u)_{j} = \alpha_{j}(\Delta u)_{j}
\end{equation*}
To evaluate the effect of the gradient limiter on the solution, discretization error consistency and order of accuracy, we decided to verify the code with gradient limiter enabled and disabled (this applies only to the second order of accuracy tests).
\end{enumerate}


